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INTRODUCTION 


Multimode  fibers  are  used  in  short-haul,  medium  bit-rate 
applications  where  their  large  core  diameters  and  numerical 
apertures  permit  the  use  of  light-emitting  diodes  as  sources  to 
realize  cost-effective  optical  systems.  These  fibers  can  guide 
several  hundred  modes  and  achieve  bandwidth-length  products  of  a 
few  gigahertz-kilometer  (references  1,  2,  3,  and  4). 

Because  fiber-optical  systems  are  not  perfectly  straight, 
uniformly  cylindrical  waveguides,  optical  power  launched  in  a 
particular  propagating  mode  does  not  necessarily  remain  in  that 
mode.  Not  only  is  the  optical  power  scattered  into  nonpropagating 
modes  and  ultimately  lost  to  the  system,  but  it  is  scattered 
between  propagating  modes  as  it  flows  through  the  components  of 
the  system.  The  characterization  of  this  power  flow  is  greatly 
complicated  by  the  fact  that  different  modes  experience  different 
group  propagation  delay  and  attenuation  through  each  segment  of 
the  system.  Consequently,  details  of  the  modification  of  the 
distribution  of  optical  power  among  the  propagating  modes  by  a 
component  are  important  in  determining  the  optical  performance  of 
any  subsequent  fiber-optical  component  in  the  system.  This 
dependency  is  known  as  the  concatenation  effect  and  should  be 
considered  in  the  specification,  design,  and  maintenance  of 
multimode-fiber-optical  systems. 
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In  this  report,  a  transfer-matrix  formalism  for  expressing  the 
effect  of  a  component  on  the  mode-power  distribution  is  discussed. 
It  is  expected  that  this  formalism  will  serve  to  both  clarify  and 
extend  matrix  approaches  to  the  description  of  this  effect  that 
have  been  presented  by  others  and  will  ultimately  result  in 
component  loss  characterization  methods  that  are  independent  of 
the  launch  conditions  chosen  to  make  the  characterization 
(references  5,  6,  7,  8,  9,  10,  11,  12,  and  13). 

BASIC  CONCEPTS 

The  mode-power  distribution  vector  u  specifies  the  power 
propagating  in  each  mode  and  is  designated  by  an  azimuthal  and 
radial  mode-number  pair.  Consequently,  a  component  can  be 
characterized  by  a  matrix  of  scattering  coefficients  s^j ,  kl  giving 
the  scattering  of  power  in  each  possible  mode  uk^  at  the  input 
into  every  possible  mode  u'^j  at  the  output.  For  simplicity,  it 
is  assumed  that  all  modes  propagate  in  the  same  direction  in  the 
system,  that  is,  reflected  power  is  considered  lost.  The  effect 
of  a  component  on  the  power  flow  can  then  be  written  as  the  matrix 
product 

K  L 

U  1 #  j  =k  l  1  x  l  1si j , kl  ukl  (!) 

where  i  =  1,...,I  and  j  =  l,...,j,  in  which  there  are  I  x  J 
possible  modes  at  the  output  and  K  x  L  possible  modes  at  the 
input . 
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Because  there  are  so  many  possible  modes,  the  scattering 
matrix,  s,  cannot  be  practically  used  to  characterize  a  component. 
Instead,  a  reduced-dimension-transfer  matrix  must  be  formed  from 
s.  Let  T  be  a  transformation  from  the  K  x  L  valued  mode-power 
distribution,  to  an  M  valued  reduced-dimension  mode-power 
distribution.  That  is,  define 


K  L 

Pm  =  X  X  ^m,kl  ukl 

K  1  1  1 


(2; 


where  m  =  1 , 


,  M,  and  p  similarly  in  terms  of  T'  and  u'  so  that 
P  Ki  =  X  Pm  ( ^ ) 


M 

X  Mnm  Pm 
m  =  1 


in  which, 


M 


nm 


=  X 


u 

X 


K 

X 


T  n,ij  sij,kl  T+ki,m 


(4) 


i  =  1  j  =  1  k  =  1  1  =  1 

defines  the  reduced-dimension  mode-power-transfer  matrix.  In  this 
definition,  the  transformation  is  the  inverse  of  Tm  in 

some  best-possible  sense,  that  is,  so  that  p'  obtained  from 
equation  3  approximates  T'u'  if  p  is  given  by  equation  2.  Here, 
the  Moore-Penrose  generalized  inverse  (see  appendix  A)  will  be 
assumed. 


Consider,  as  an  example  of  dimension  reduction,  the 
transformation  from  the  full  radial/azimuthal  mode  description  to 
the  reduced-dimension  fundamental-mode-number  description  in  which 
the  LP  modes  are  grouped  according  to  a  fundamental  mode  number 
m  =  21-1+k  (references  1,  14,  and  15).  While  the  actual  form  of 
the  transformation  needed  to  accomplish  this  reduction  in 
dimension  is  not  necessarily  specified,  evidently  there  is  very 
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little  error  incurred  in  its  use.  Therefore,  it  must  hold,  in  the 
appropriate  sense  (appendix  A) ,  that 

.  £  .  Z  Tn,ij  T  ij,m  =  ^nm 

1  =  1  j  =  i  J 

as  well  as 

Z  T+ij/ltl  TM(kl  =  5ij  Sn  (6) 

m  =  1 

where  is  the  Kronecker  delta.  In  what  follows,  it  will  be 

assumed  that  the  fundamental  mode  groups  of  the  linearly  polarized 
modes  form  the  basis  for  further  reduction  in  dimension.  This  is 
appropriate  for  loss  characterization  because  while  these  modes 
may  not  be  strongly  coupled  within  a  group  (references  16  and  17), 
they  appear  to  be  essentially  monotonic  with  respect  to  loss,  that 
is,  higher  order  modes  are  generally  more  attenuated  than  lower 
order  modes. 

MODE-BLOCK  REDUCTION 

An  explicit  example  of  a  transformation,  T,  that  combines 
fundamental-mode  groups  into  mode  blocks  is  the  matrix 

1  1...1  0  0...0  0  0 ...  0 ...  0  0...0 

0  0...0  1  1...1  0  0...0...0  0...0  (7) 

0  0...0  0  0...0  U  0 ...  0 ...  1  1...1 

where  there  are  M  rows,  one  for  each  block  and  N  columns,  one  for 
each  fundamental-mode  group,  and,  if  the  partition  of  the  mode 
groups  is  Nm  in  the  mth  mode  block,  the  number  of  ones  in  the 
diagonal  non-zero  partition  of  the  mth  row  of  T  will  be  Nm.  Here, 
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the  modes  are  all  weighted  either  zero  or  one.  By  inspection,  the 
generalized  inverse  of  this  transformation,  T+,  is 


1/N2  0  ...  0 

1/N-l  0  ...  0 

1/N-l  0  ...  0 

o  i/n2  . . .  o 

0  i/n2  ...  o 

o  i/n2  ...  0 

0  0  ...  0 

0  0  ...  0 

0  0  ...  0 

0  0  ...  i/nm 

0  0  i/nm 


I  0  0  ...  1/Nm  I 

where  there  are  N  rows  and  M  columns,  in  which  there  are  Nj„ 
non-zero  elements  in  the  diagonal  partition  mth  column. 


(8) 


The  product  TT+  is  the  M  by  M  identity  matrix 

1  0  ...  0 
0  1  ...  0 

00  ...  1 

while  the  product  T+T  is  the  N  by  N  matrix 


Ei/Ni  0  ...  0 

o  e2/n2  ...  0 

0  0  .  .  .  Ejij/Njij 

in  which  Em  is  the  Nm  by  Nm  matrix  of  all  ones 


(10) 


1  1  . . .  1 

1  1  . . .  1 

1  1  .  . .  1 


(ID 


Clearly,  T+  satisfies  the  conditions  to  be  the  unique  generalized 
inverse  of  T  and  is,  in  the  Frobinius  norm,  the  matrix  that  makes 
T+T  closest  to  the  N  by  N  identity  matrix  (see  appendix  A) . 
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A  generalization  of  mode-block  reduction  is  to  choose  T  so 

that 


N 

I  wn  Tin  Tjn  =  5ij 
n  =  1 

with  i , j  =1 , . . . , M .  Consequently , 

T+kl  =  W1  Tlk 

with  l,k=l,...,N.  In  this  instance,  TT+  is  again  the  M  by  M 
identity  matrix,  and  the  N  by  N  matrix  T+T  with  elements 

M 

(T+T) ij  =  I  wm  Tmi  Tmj 

is  in  general  not  the  identity  matrix  unless  M  =  N.  This 
extension  amounts  to  taking  the  rows  of  T  to  be  the  first  M 
polynomials  from  an  N  dimensional  set  of  discrete  orthogonal 
polynomials  having  a  weight  function  wn. 


These  examples  illustrate  the  result  that  when  the 
reduced-dimension  transfer  matrix  for  a  concatenation  of  two 
components  is  expressed  in  terms  of  the  individual  components  in 
the  form 

M12  =  M 2  M1  (15) 

which,  by  equation  4,  can  be  rewritten  as 

m12  _  Tns2  t-+  T'  S1  T+  (16) 

it  is  not  identical  to  the  reduced-dimension  transfer  matrix  for 
the  same  system  that  would  be  formed  directly  from  the  product  of 
the  two  scattering  matrices  s2s1  because  of  the  presence  of  the 
product  T'+T  ,  which,  in  general,  is  not  an  identity  matrix. 
Consequently,  some  error  must  always  be  expected  if  a  system  is 
described  in  terms  of  the  reduced-dimension  transfer  matrices  of 
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its  components;  however,  the  error  will  be  minimal  in  the 
Frobinius  norm  if  T+  is  taken  to  be  the  Moore-Penrose  generalized 
inverse  of  T  (see  appendix  A) . 


CONTINUOUS  MODE-POWER  DISTRIBUTION 


Measurement  considerations  prohibit  the  determination  of  even 
the  fundamental-mode-number  distribution  so  that  further 
simplification  of  the  mode  distribution  must  be  introduced.  If 
the  propagating  field  intensity  is  axially  symmetric  and  the 
spectral  width  of  the  source  51  satisfies  a  minimum  width 
requirement,  which  for  a  parabolic  index  profile  is 


X  a  k  ngw 


where  X  is  the  free-space  wavelength,  k  is  the  free-space  wave 
number,  A  is  the  relative-refractive-index  difference,  a  is  the 
core  radius,  and  n^  is  the  effective-group  index  of  refraction, 
there  will  be  many  modes  excited  and  they  will  be  spaced  more 
closely  than  their  spectral  width.  In  this  situation,  the 
mode-power  distribution  can  be  taken  to  be  a  continuous  function 
of  a  continuous  mode  parameter  rather  than  a  discrete  function  of 
a  discrete  mode  index  (references  15  and  16) .  More  important 
here,  there  exists  a  relationship  between  the  continuous 
approximation  to  the  mode-power  distribution  and  the  measurable 
near-field  pattern,  I  (p) ,  on  the  end  face  of  a  fiber  (references 
16,  18,  19,  20,  21,  22,  23,  and  24).  This  relationship  takes 
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different  forms  depending  on  the  choice  of  mode  parameter,  so  some 
consideration  of  the  mode  parameter  is  necessary. 


In  the  continuous  approximation,  it  can  be  written  that  the 
total  power  propagating  W  is 


N 

W  =  £  pn  =  N 

n  =  1 


pz(z)dz 


(18) 


where  pn  is  the  optical  power  in  the  nth  fundamental-mode  group,  N 
is  the  number  of  such  mode  groups,  and  pz(z)  is  the  continuous 
approximation  to  the  discrete  distribution  pn.  Here,  z  is  the 
continuous  extension  of  the  discrete  parameter  n  normalized  so 
that  0  <  z  <  1  and  chosen  so  that  the  modes  have  a  uniform  density 
distribution  N.  Other  choices  for  the  continuous  mode  parameter, 
for  example,  x:z=g(x),  where  g(0)  =  0  and  g(l)  =  1  so  that  0  <x< 

1,  yield  the  generalization 


w  = 


ax<x)  Px(x>  dx 


(19) 


where 


and 


Px(x)  =  Pz[<3(x)  3 


(20) 


in  which 


g(x) 


1 

N 


'X 

ax (x) dx 

0 


(21) 


N 


*1 

ox (x) dx 

v  0 


(22) 


We  shall  call  px(x)  the  ax (x) -density  distribution  corresponding 
to  the  continuous  mode  parameter  x=g-1(z).  These  transformed 
descriptions  of  the  mode  distribution  are  all  equivalent; 
furthermore,  even  though  the  uniform  density  distribution  seems  to 


8 


have  a  more  natural  physical  interpretation  in  terms  of  the 
discrete  mode  distribution,  it  is  not  the  most  useful  because  the 
connection  between  a  continuous  mode-power  distribution  and  the 
emittance,  I  (p),  on  the  face  of  the  fiber  has  not  been  cast  in 
terms  of  the  uniform-density  distribution. 

The  relationship  between  I (p)  and  px(x)  shall  be  written  here, 
neglecting  leaky  modes,  as 

v2  r1 

I(S)  =  - -  px(x)dx  ,  s  =  p/a  (23) 

2jra  J  f  ( s) 

where  V  is  the  normalized  frequency,  f(s)  is  the  index  of 
refraction  profile  tacitly  assumed  here  to  be  monotonic,  and  p  is 
the  radial  coordinate  at  the  end  face  of  the  fiber.  The 
normalized  frequency  can  be  written 

V  =  a  k  n(0)  /2a'  (24) 

where  n(0)  is  index  of  refraction  at  the  center  of  the  fiber  core 
and  the  index  of  refraction  profile  is  defined  by 

n2(s)  =  n2 ( 0)  [1  -  2A f (s )  ]  (25) 

and  the  condition  that  0  <  f(s)  <1. 

By  differentiating  both  sides  of  equation  23,  an  explicit 
expression  for  px(x)  is  obtained 

V2 

l  (s)  = - r  px|_f  (s)Jf  *  (s)  .  (26) 

2na* 

This  expression  is  the  one  generally  used  to  determine  the 
mode-power  distribution  from  measurements  of  the  near-field 
pattern  (references  8,  18,  23,  and  24).  Because  equation  26 
involves  the  ratio  of  the  derivatives  of  the  two  quantities  I(s) 
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and  f(s)  to  determine  px(x),  the  presence  of  noise  or  distortion 
in  the  measurements  of  these  quantities  can  introduce  a  large 
uncertainty  into  px(x)  particularly  near  extrema  in  f(s).  An 
alternative  method  using  equation  23  directly  will  be  presented 
later  in  the  section  on  experimental  determination. 

Another  frequently  occurring  form  for  the  relationship 
between  the  I (p)  and  the  mode  distribution  uses  the  parameter 
y=Vx  (references  12,  16,  and  25).  For  this  change  of  variable, 
equation  23  yields 

fl 

I(S)  =  — -  y  Py(y)  dy  (27) 

n  a2 

yfW 

and  equation  26  takes  the  form 

V2  _ 

I'(s)  - -  py  [yfT¥T]  f'(s)  .  (28) 

27ra2  Y 


That  neither  x  nor  y  is  the  uniform  density  parameter  z  can  be 
seen  from  the  total  power  flow  in  the  fiber,  which  can  be  written 
in  the  case  of  x,  for  example,  either  as 


W  =  2na 


s 


I(s) 


ds 


(29) 


or,  by  introducing  equation  23  and  changing  the  order  of  the  s  and 
x  integration,  as 


W 


v2  r1 

-  [f-1(x) ]2  px  dx 

2  Jo 


(30) 


from  which  it  can  be  seen  that  the  mode  density  corresponding  to 
the  parameter  x  is 
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ax(S) 


V2 

2 


[f1(x)  ]2 


(31) 


A  more  general  example  that  is  useful  for  computational 
purposes  is  that  of  the  power  profile  where  f(s)  =  sa  and  g(y)  = 
yP .  In  this  case, 

bv2  r1 

I(S)  =  - -  y^-1  Py(y)  dy  (32) 

27ra2  Jsa//3  y 

and 

CTy(y)  =  - -  yP  ( 3-+2/0! )  ~1  .  (33) 


Now  it  is  seen  that  if  p  is  chosen  to  have  the  value  (3  =  a/ (a+2)  , 


then  we  have  the  uniform-density  choice  of  the  mode  parameter, 
that  is,  y  =  z  with  the  uniform  mode  density 


V2  a 
2 (a+2 ) 


(34) 


being  the  number  of  propagating  modes  N.  The  near-field  pattern 
is  then  given  by 

v2 

I(s)  =  - — -  z-2/(a+2)  pz(z)  dz  .  (35) 

2?ra2(a+2)  J  sa+2 

Consequently,  we  obtain  for  the  uniform-density  mode-power 
distribution 


Pz(z)  =  " 


Z  (!”<*)/  (2+Q)  J  '  (g) 


s  =  Z ( a+2 ) 


•  (36) 


The  factor  that  appears  in  front  of  the  integral  in  equation 
23  may  take  different  forms  depending  on  the  choice  of  units  of 
px(x) ,  which  need  not  be  power  per  mode  as  we  have  selected.  See 
equation  19,  where  ox(x)  has  the  units  of  modes  per  unit 
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continuum-mode  parameter.  For  example,  choosing  units  of 
radiance,  that  is,  W/cm^sr  places  quite  a  different  physical 
interpretation  on  px(x) ,  but  nevertheless  results  in  a  density 
function  ax(x)  having  a  squared  dependence  on  the  inverse  profile 
function  just  as  in  equation  31  (references  1  and  20) .  As  will 
become  evident  later,  the  essential  feature  of  equation  23  is  that 
it  is  a  linear  relationship  between  I(s)  and  px(x). 


MATRIX  FORMALISM  FOR  THE  CONTINUOUS  MODE  DISTRIBUTION 


Just  as  in  the  discrete  case,  a  matrix  formalism  related  to 
the  measurable  near-field  patterns  can  now  be  introduced  by 
assuming  there  exists  a  linear  relationship  of  the  form 


P'z<z) 


'1 

Pz(2'>hz(z' >z)  dz' 

o 


(37) 


where  the  primed  and  unprimed  p  designate,  respectively,  the 
output  and  input  mode-power  distributions  of  the  element,  and  the 
kernel  h(z',z)  is  a  system  response  function  that  gives, 
essentially,  the  fraction  of  the  power  propagating  in  the  region 
dz 1  about  z'  scattered  into  the  region  dz  about  z. 


A  matrix  representation  of  h(z’,z)  is  obtained  by  expanding 
p'z(z)  and  pz(z')  in  a  suitable  set  of  orthogonal  functions 
{<*»n(z))'  where 

>1 

w(z)  0i(z)  0j (z)  dz  =  6 ± j  (38) 

*  0 

is  the  Kronecker  delta,  and  w(z)  is  the  weight  function  of  the 
set  ( <^m ( z )  ) .  Thus, 
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co 


Pz(z') 

=  I  cm  <*>m(z') 

m  =  l 

(39) 

00 

Pz(z) 

=  l  c'm  ^m(z)  • 

m  =  1 

(40) 

Substituting  equations  39  and  40  into  equation  37,  multiplying 
both  sides  of  the  result  by  w(z)  0r(z),  and  integrating  over  z 
yields  upon  application  of  equation  38 


00 


I  Hr,mcm'  r  ”  i/2» • • • 
m  =  1 


(41) 


where 


1 

*1 

" 

•1 

a 

n 

3 

ii 

w(z)<*>r(z) 

0m(z')h(z',z)dz' 

J 

0 

_  \ 

0 

(42) 


are  the  elements  of  the  desired  matrix  representation  of  h(z',z). 
Practically,  the  summations  in  equations  39  and  40  will  have  to  be 
truncated  at  some  term  N  with  the  results  that  H  will  be  an  N  by  N 
matrix  approximation  to  h(z',z)  and  c'  and  c  will  be  N  by  1  matrix 
(vector)  approximations  to  p'z(z)  and  pz(z),  respectively.  In 
matrix  notation,  equation  41  becomes 

c'  =  H  c  .  (43) 


While  the  optical  source  can  be  represented  by  the  vector  c, 
the  representation  of  the  photodetector  requires  the  additional 
assumption  that  the  signal  current  j  can  be  obtained  by  a  simple 
weighting  of  the  mode-power  distribution  at  the  input  of  the 
photodetector.  Thus, 


j 


pz(z)  Rz(z)  dz 


(44) 


where ,  now 
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W(z)  00 

Rz(z)  =  -  l  rm  ^m(z)  (45) 

N  m  =  1 

yielding  upon  application  of  equation  38 

3  •  X  rm  cm  (46> 

m  =  1 

where  the  rm  are  the  vector  coefficients  of  the  photodetector 
response  function  Rz(z).  This  expansion  must  be  truncated,  so  r 
becomes  a  1  by  N  matrix  (vector)  and  the  scalar  photocurrent  j  can 
be  approximated  as  the  dot  product 

j  =  r  •  c  .  (47) 

Using  a  relationship  such  as  equation  43,  equation  47  can  be 
rewritten  as 

j  =  r  •  H  c  .  (48) 

This  is  the  basic  system  equation  in  the  transfer-matrix  formalism 
and  accounts  for  the  mode-power  distribution  of  the  source,  c,  the 
intermode  scattering  in  the  optical  system,  H,  and  the  selective 
response  of  the  photodetector,  r. 

The  details  of  the  form  chosen  for  equation  44  are  unimportant 
for  our  purposes,  and  an  equation  like  equation  47  could  have  been 
derived  from  consideration  of  the  discrete  formalism  as  well. 

EXPERIMENTAL  DETERMINATION 

Neither  equation  4  nor  42  is  of  much  use  other  than  to  provide 
a  theoretical  background  because  neither  the  matrix  s^j  nor  the 
continuum  function  h(z',z)  appearing  in  them  are  likely  to  be 
obtainable.  An  exception  to  this  might  occur  in  situations  where 
classes  of  components,  such  as  lengths  of  fiber,  can  be 
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parametrized  with  only  a  few  parameters  resulting  in  significant 
simplification  in  the  required  measurements  (references  26  and 
27) .  Therefore,  H  must  be  determined  directly  from  a  set  of  input 
and  output  radiometric  measurements  from  which  input  and  output 
distributions  px(x)  can  be  found.  In  the  mode-continuum 
approximation,  equation  26  provides  such  a  needed  relationship 
between  the  mode  distribution  and  measurable  near-field  patterns. 
To  avoid  the  numerical  problems  involved  with  evaluating 
derivatives  of  noisy  data,  as  well  as  possible  distortion-induced 
singularities,  the  integral  form,  equation  23,  can  be  used  in  a 
least-square  approach  as  follows. 


Begin  by  expanding  the  mode-power  distribution  in  a  finite 
series 


P 

Px(x)  =  I  cm  <*>m(x>  (49) 

m  =  1 

where  P  is  small  compared  to  the  number  of  propagating  modes,  but 
must  be  equal  to  or  larger  than  N,  the  dimension  of  the  desired 
transfer  matrix,  and  {<£m(x)  )  is  a  suitable  set  of  orthogonal 
functions.  Now,  in  addition  to  errors  produced  by  noise  in  the 
measured  quantities,  there  will  be  an  error  caused  by  the 
truncation  of  the  expansion  at  N.  It  is  possible  that  the 
truncation  error  could  be  reduced  if  a  rational  function  were 
used  rather  than  a  series  expansion;  however,  here  we  are  not 
interested  in  merely  finding  the  mode-power  distribution  px(x), 
but  in  obtaining  a  vector  representation  of  it  (refernce  28). 
Therefore,  in  order  to  reduce  the  effects  of  all  errors,  the  cn 
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are  chosen  to  minimize  the  mean-square  error  in  the  power  flow 
defined  by 


^  =  2na‘ 


1  r 


I(s)  - 


V2 


l  c: 


m 


0m(x)dx 


f  (s) 


2ma2  m  =  1 

The  normal  equations  for  the  cn  are  obtained  from 

3  E2 


ds 


(50) 


=  0 


3c, 


r  =  1,  2, 


(51) 


which  can  be  written  in  matrix  form 


A  c  =  b 


(52) 


where 


rl 


s  i  ( s)  0j(s)  ds  i ,  j  =  1,  ...,p  (53) 


s  I(s)  0 j (s)  ds  j  =  1, 


(54) 


and 


rl 


’/’n  (s)  = 


2ira‘ 


f  (s) 


0n(x)dx  ,  n  =  1,  . . . ,p 


(55) 


The  cn  obtained  from  the  solution  of  the  normal  equations  can 
be  used  directly  to  represent  the  reduced  dimension  mode-power 
distribution  simply  by  truncating  the  vectors  after  N  components, 
or  the  further  reduction  in  dimension  can  be  obtained  by  a 
subsequent  transformation.  For  example,  to  obtain  the  N 
components  pj  for  a  mode-block  representation,  the  following 
transformation  could  be  applied  in  the  case  where  P  is 
significantly  greater  than  N 

P 

pj  =  l  Tjn  cn  j  =  1,  . . . ,  N  (56) 

n  =  1 

where 
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and 


:n 


“j  = 


x-i_ 


ax(x)  <*>n(x)dx,  x0  =  0,  xn  =  1 


D-l 


Xj.i 


ax(x)  dx/ 


ax(x)dx  j  =  1,  ...»  N 


(57) 


(58) 


define  (Xj)  in  terms  of  the  mode-block  partition  {m-^}  so  that  the 
resulting  transfer  matrices  would  be  independent  of  the  choice  of 
mode  parameter  x. 


If  a  set  of  N  measurements  of  input  and  output  near-field 
patterns  are  made  for  an  element  of  a  system,  then  the  N  solutions 
to  equation  52  for  the  input  and  output  mode-power-distribution 
coefficients,  or  some  transformation  of  them,  can  be  used  in 
equation  43  to  find  H,  the  transfer  matrix  for  that  element. 
Introducing  the  notation  for  an  N  by  N  matrix 

V  (a)  =  (v-l  |  v2  I  •  .  •  |vN)  (°)  (59) 

where  are  N  by  1  column  vectors  and  a  =  0,  or  1  indicates  input 
or  output,  respectively,  the  formal  solution  for  H  can  be  written 
H  =  [a(1> j-1  (b-J  . . . |bN) U>  [(b1|...|bN)  (°) ] “I  A ( °)  (60) 

which  requires,  not  only  that  A  have  an  inverse,  but  that  the 
matrix  made  up  of  the  vectors  bj^  obtained  by  equation  54  from  the 
N  input  near-field  patterns  have  an  inverse,  that  is,  that  these  N 
vectors  be  linearly  independent.  This  requirement  to  have  N 
independent  input  conditions  places  the  practical  upper  limit  on  N 
and,  consequently,  the  accuracy  of  reduced-dimension  transfer 
matrix  methods. 


Similarly,  N  independent  measurements  of  the  response  of  a 
photodetector  can  be  used  to  form  a  matrix  equation  whose  formal 
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(61) 


solution  for  the  response  vector  is 

r  =  J  [(b1|...|bn)  (0)]-l  A(0) 

where  J  is  the  row  vector  formed  from  the  N  photocurrent 
measurements . 

While  these  equations  represent  solutions  for  H  and  r,  there 
are  advantages  obtained  by  not  using  them  for  numerical 
computation.  Rather,  a  more  robust  technique  such  as  a  singular 
value  decomposition  of  the  input  data  matrix  made  up  of  the 
near-field  pattern  vectors  should  be  used.  Thus,  rewriting  the 
data  matrix  in  terms  of  its  singular-value  decomposition 

D(°)  =  (Cil ... |cK) (°)  =  USVT  (62) 

where  now,  K,  the  number  of  measurements,  may  be  greater  than  N, 
the  dimension  of  the  desired  reduced  dimension  transfer  matrix. 

Not  all  K  measured  vectors  need  be  independent.  Now,  from  the 
pseudoreverse  VS+UT  of  USVT,  a  least-square  fit  to  H  can  be 
found 

H  =  (CjJ  . .  .  |cK)  (!)  VS+UT  (63) 

without  calculating  explicitly  [D^0)  d(°)t]~1.  Furthermore,  from 
consideration  of  the  singular  values  contained  in  the  diagonal 
matrix  S  a  condition  number  for  d(°)  can  be  obtained  that  can  be 
used  to  evaluate  the  magnitude  of  computational  and  random  errors 
that  might  occur  in  H  as  a  result  of  round  off  error  in  the 
computation  and  noise  in  the  input  and  output  data  matrices.  More 
importantly,  the  condition  numbers  can  be  used  to  ascertain  the 
degree  of  independence  of  at  least  N  of  the  K  measured  input  data 
vectors. 
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Similar  considerations  militate  towards  the  singular  value 
decomposition  of  the  input  data  matrix  to  obtain  a  least-square 
fit  to  the  detector  response  vector  from  the  output-current  vector 
J.  In  both  cases,  the  input  data  matrices  made  up  of  the  c  vector 
columns,  should  be  obtained  by  some  robust  technique,  such  as 
Cholesky  decomposition  of  the  symmetric  A  matrix  applied  to  the 
primary  data  vector  b,  rather  than  direct  inversion  of  A. 

MEASUREMENT  CONSIDERATIONS 


The  general  conditions  under  which  the  measurements  should  be 
made  is  that  the  apparatus  remain  optically  and  electronically 
stable  for  a  given  launch  configuration  during  the  time  it  takes 
to  measure  a  pair  of  input  and  output  near-field  patterns  for  a 
component,  the  output  current  and  input  near-field  pattern  for  a 
detector,  or  a  sequence  of  output  near-field  patterns  of  a  source. 
The  near-field  patterns,  I(sa),  are  the  radiant  emittance  on  a 
fiber  end  face  and  they  have  the  units  of  power  per  unit  area. 
Typically,  near-field  patterns  are  determined  from  a  relative 
measurement,  M(u) ,  of  the  radiant  emittance  of  an  enlarged  image 
of  the  fiber  end  face,  an  absolute  measurement  of  the  total 
radiant  flux  from  the  fiber  end,  W,  and  a  measurement  of  the  fiber 
core  area,  ir a2.  Thus, 

I (as)  =  M(bs)/G  (64) 

here  the  instrument  gain,  G,  is  defined  to  include  the  image 
magnification  b/a,  where  b  is  the  radius  of  the  core  image,  as 
well  as  the  electronic  amplification.  The  gain  is  given  by 
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G 


(65) 


jra 


2 


W 


•b 

u  M(u)du 

J  o _ 

•b 

u  du 
o 


Usually,  the  image  coordinate,  u,  is  scaled  so  as  to  make  b 
approximately  full  scale. 


The  input  radiometric  quantity  used  to  characterize  optical 
detectors  can  be  either  flux  (power)  or  irradiance  (power  per  unit 
area) ,  depending  on  the  intended  use  of  the  detector.  If  the 
detector  is  to  be  irradiated  with  a  spot  of  radiation  that  does 
not  fill  the  detector,  flux  is  generally  chosen.  If,  on  the  other 
hand,  the  detector  is  to  be  overfilled  with  relatively  uniform 
radiation,  irradiance  would  be  used.  For  fiber-optic  systems,  the 
detector  can  be  expected  to  be  underfilled,  so  that  flux  would  be 
the  expected  choice.  However,  most  measurements  are  related  to 
the  determination  of  transfer  matrices  of  components,  and  if 
irradiance  is  chosen  to  be  the  quantity  characterizing  the  system, 
it  is  not  necessary  to  have  values  for  either  the  first  factor  of 
G,  equation  65,  or  the  normalized  frequency  V,  equation  24,  to 
determine  transfer  matrices.  This  is  because  choosing  irradiance 
rather  than  flux  to  characterize  the  power  flow  in  the  system  is 
equivalent  to  using  V2p(x)/27ra2  rather  than  V2p(x)/?  to 
characterize  the  mode-power  distribution.  In  any  case,  if  the 
instrument  gain  were  to  be  changed  between  the  input  and  output 
near-field  measurements  of  a  pair,  the  gain  ratio  of  the  two 
measurements  would  have  to  be  determined. 
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Of  course,  whatever  the  choice  of  characterizing  radiometric 
quantity,  in  order  to  determine  source  or  detector  response 
vectors,  the  gain,  G,  is  needed  requiring  measurements  of  both  the 
core  radius,  a  and  the  power,  W.  If  an  absolute  determination  of 
the  mode-power  distribution  is  desired,  then,  in  addition,  a  value 
of  the  V  parameter  must  be  available. 

The  transfer  matrices,  as  well  as  the  source  and  detector 
response  vectors,  will  be  different  by  scalar  factors  for 
different  choices  of  normalization  units.  For  example,  in  the 
mode-block  case,  matrices  determined  with  irradiance  as  the 
significant  quantity  must  be  multiplied  by  the  ratio  of  the  areas 
of  the  output  to  the  input  fiber  cores  in  order  to  obtain  the 
matrices  corresponding  to  flux  as  the  significant  quantity. 
Consequently,  for  irradiance  based  transfer  matrices,  the  upper 
limits  for  matrix-element  and  column-sum  size  become  the  ratio  of 
input  to  output  fiber  core  areas  rather  than  one  as  it  is  for  flux 
based  transfer  matrices  (see  TRANSFER  MATRIX  INTERPRETATION) . 

ALTERNATIVE  VIEWPOINT 

In  addition  to  the  mode-block  and  expansion-coefficient 
formalisms,  other  approaches  to  generating  sets  of  consistent 
transfer  matrices  are  possible.  For  example,  the  matrix  A  defined 
in  equation  53  can  be  considered  to  be  a  transformation  of  the 
vector  c  into  the  vector  b  defined  in  equation  54.  In  such  a 
representation,  the  measured  transfer  matrix  becomes 
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(66) 


H  =  (bjJ  ... | (bK) (1)  t (bxl  —  |bK) (0)1 
and  the  measured  detector-response  vector  becomes 

5  =  5  [<5,|...|bK)«»]  +  .  (67) 

In  this  representation,  neither  the  matrix  A  nor  its  inverse  needs 
to  be  calculated. 


Even  greater  simplification  results  if  a  transformation  T  on  b 
can  be  found  such  that 


<5i  “  .  I  xTij  bj 

for  which 


(68) 


N 

£i(s)  =.  I  Tij  *j(s) 

1  *  1 


(69) 


so  that  from  equation  54 


Si  = 


sl(s)  ei(s)  ds 


(70) 


where  (£i(s)}  is  a  set  of  orthogonal  functions  on  the  domain 
(0,1)  with  weighting  function  s.  In  this  case,  the  set  (q^)  is 
the  expansion  coefficients  of  I(s)  for  the  chosen  set  of 
orthogonal  functions,  for  example,  Jacobi  polynomials,  and  are 
used  in  place  of  the  (b^)  to  calculate  transfer  matrices.  The 
transformation  T  and  functions  0 j_  are  not  explicitly  required; 
however,  their  existence  must  be  assumed.  It  has  not  yet  been 
verified  either  theoretically  or  experimentally,  but  it  seems 
plausible  that  they  do  exist  from  consideration  of  the 
differential  form  of  equation  55  and  the  idea  that  sets  of 
orthogonal  functions  can  be  found  whose  derivatives  are  also  sets 
of  orthogonal  functions. 
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An  advantage  of  this  extension  is  that  the  index  of  refraction 
function,  f(s),  does  not  necessarily  enter  the  calculation  unless 


the  mode-power  distribution,  p(x) ,  is  required,  and  the  core 
radius,  a,  enters  only  relatively  in  the  definition,  s  =  p/a . 
Furthermore,  if  the  {£j.(s)>  are  chosen  to  be  rectangular  functions 
that  are  one  for  0  =  s0<s1 . . . <s^_^<s<s^<s^+1 . . . <sN  =  1  and  zero 
otherwise,  then  the  are  directly  measurable  as  the  power  within 
an  annulus  in  the  near-field  pattern.  An  appropriate 
choice  for  the  partition  {s^>  might  be  in  this  case 


s[ l-f (s) ]ds  =  — 
Si_i  N 


s  [  l-f ( s) ] ds 


reintroducing  f(s)  into  the  calculation  of  q. 


(71) 


EQUIVALENCE  OF  REPRESENTATIONS 


The  various  representations  of  the  reduced-dimension  transfer 
matrices  can  be  related  to  one  another  through  the  full  matrix  of 
scattering  coefficients,  s,  and  the  dimension-reducing 
transformations ,  T,  corresponding  to  each  of  the  reduced-dimension 
representations.  Let 

H  =  TsT+  (72) 

be  an  M-dimensional  representation  of  s  and 

ft  =  (73) 

be  an  M-dimensional  representation  of  s,  then,  in  a  best-possible 
sense , 

H  a  T'fr+ftfa’+  .  (74) 
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This  relationship  is  only  approximate  because,  while  H  and  H, 
respectively,  are  defined  by  equations  72  and  73,  it  is 
generally  true  that 

S  x  fr+frfc  (75) 

because,  in  general,  the  information  lost  in  conversion  of  s  into 

ft  cannot  be  fully  recovered.  Consequently,  equation  74  should  be 

,  -v 

most  useful  if  MsM. 


For  example,  consider  the  transformation  from  a  mode-block 
representation  with  ft  =  3  to  one  which  has  M  =  2,  where  the 
partitions  are  related  by  the  restriction  m1<mj<m^+m2-  In  this 
instance,  equations  7  and  8  yield 


(76) 


and 


(77) 
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so  that,  according  to  equation  74, 


H 


*  (ft. 


11  ~  \nll 
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A/  A/  A/  A/ 

m2-m3  \  \  ml  m2-In3  v  <v  mj^-mj^ 
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A/ 
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(81) 


In  contrast,  the  generalized  mode-block  reduction  defined  by 

.  -v 

equations  2,  12,  and  13  provides  in  the  case  for  which  MsM  and  the 

same  set  of  discrete  polynomials  is  used  for  both  representations 
that  Hjj  ~  fti j ;  that  is,  H  is  a  principal  diagonal  matrix  of  ft. 
This  simple  relationship  is  a  result  of  the  orthogonality  of  the 
basis  functions  for  the  two  representations.  Furthermore,  by 
letting  M  be  reduced  to  one,  it  follows  that,  in  this  case, 
should  be  the  scalar  quantity  best  characterizing  the  transmission 
of  power,  or  power  per  unit  area  as  the  case  may  be,  through  a 
component.  In  the  mode-block  case,  the  analogous  scalar 
transmission  for  a  component  is  given  by  the  transformation 

T^+  =  ( 1 ,  !,...!)  (82) 
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(83) 


*T+  =  (-, 

A,  A/ 

M  M 


A/  t  A/  # 

of  the  M  dimensional  matrix  H  to  the  one  dimensional  matrix  Hu- 
This  results  in 


(84) 


In  general,  the  matrix  T  that  transforms  s  into  H  will  not  be 
available,  so  an  alternative  approximate  transformation  between 
representations  is  desirable.  Suppose  that  the  power-flow  vectors 
for  the  two  representations  can  be  related  by  a  transformation  Q 
such  that  c  =  Qc.  If  c '  =  He  and  c '  =  ftc,  then 

Q+Q  c'  «  Q+Q  H  Q+  QC  (85) 

or 

C'  »  QH  Q+  c  (86) 

where  Q+  is  the  generalized  inverse  of  Q,  defined  so  that  c  w  Q+c 
so  that  Q+Q  «  I.  Hence, 

H  «  QH  Q+  (87) 

that  is,  Q  and  Q+  replace  TT+  and  TT+,  respectively.  Now  if  the 
mode-power  distribution  can  be  approximated  either  by 

M 

P (x)  *  l  Ci  0i(x)  (88) 
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(89) 


p(X)  a  l  Ci  </>i(X) 

i  =  1 


where  {<p±)  and  (0^)  are  sets  of  orthogonal  functions  such  that 


r  l 

w(x)  0i  0j  ( x)  dx  6  i  j 
J  0 

1  w { x)  ^i(x)  (x) dx  "  6i3 

J  0 


It  then  follows  that 


r1  2  a 

Ci  w(x)  0i (x) dx  a  l 
Jo  j  =  1 

so  that 


w(x)0i (x)0j (x) dx  Cj 


w(x)0i(x)0j (x) dx 


-iJ  ~  ri 


w(x)0i2 (x) dx 


Similarly , 


+  Jo 

*  71 


(x) dx 


'1 

w(x)  ^i2(x)dx 
J  0 


Transformations  within  sets  of  basis  functions,  that  is,  where 
{<^i>  =  { 0  i )  ,  yield  immediately  the  simple  result  that  ft  ^  j  a  Hi  j  , 
i , j  =1,  . . . ,  M,  where  M<ft,  just  as  previously  obtained. 


In  the  case  of  a  mode-block  representation,  the  {0i)  are  the 
rectangular  functions  and  w(x)  is  proportional  to  the  mode 
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density,  a  (x) .  The  case  of  a  transformation  between  two 
mode-block  representations  can  be  simplified  if  the 
representations  have  been  defined  in  terms  of  a  fractional  number 
of  modes  per  block,  as  in  the  use  of  equation  58,  so  that  the  mode 
density  a (x)  can  be  taken  to  be  uniform  even  though  it  may  not 
have  been  in  the  original  experimental  determination,  as  in  the 
use  of  equation  23.  In  this  situation,  the  Q^j  represent  the 
overlap  of  the  {m^}  blocks  on  the  {m^)  blocks  and  the  reverse  for 
j ,  and  the  use  of  diagrams  such  as  those  shown  figures  1  and  2 
can  aid  in  evaluating  this  overlap.  For  example,  consider  a 
transformation  from  M  =  2  to  ft  =  3 ,  where,  for  simplicity,  the 

Ay  'V 

mode-partitions  satisfy  the  restriction  m1>m1+m2.  Aided  by  figure 
1,  we  obtain  from  equation  93 
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Figure  1 .  Aid  for  obtaining  the  matrix  elements  in  equation  95. 


Figure  2.  Aid  for  obtaining  the  matrix  elements  in  equation  96. 
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here  as  well  as  in  the  earlier  example,  the  column  sums  of  the 
transformation  matrices  add  to  one. 

Although  we  have  shown  that  all  representations  are 
equivalent  in  the  sense  that  one  can  be  obtained  approximately 
from  another,  it  must  be  kept  in  mind  that  each  transformation 
causes  some  of  the  information  that  was  contained  in  the  original 
measurements  of  the  near- field  patterns  to  be  lost.  Consequently, 
whenever  it  is  possible,  original  data  should  be  used  or  augmented 
to  make  new  calculations  rather  than  attempting  to  transform 
previously  determined  matrices. 


TRANSFER  MATRIX  INTERPRETATION 


Except  in  the  mode-block  approach,  the  basis  functions  for  the 
expansion  of  the  mode-power  distribution  are  oscillating 
functions;  consequently,  there  is  no  simple  general  physical 
interpretation  of  the  vector  components  or  of  the  transfer-matrix 
elements,  which  can  be  either  positive  or  negative  and  of 
unpredictable  magnitude.  Discussion  of  the  mode-block 
interpretation  is  postponed  until  following  the  mention  of 
two  special  cases  in  which  the  lowest  order  component  of  a  state 
vector,  c,  can  be  related  to  the  total  optical  power  flow. 

The  first  case  is  obtained  by  choosing  a  representation  for 
which  c  is  proportional  to  q,  defined  by  equations  68  to  70,  and 
noting  that  the  lowest  order  function,  is  a  constant.  The 
second  case  arises  if  the  weighting  function,  w(x),  of  the  {0n(x)> 
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is  chosen  to  be  the  mode  density,  ax(x) ,  defined  by  equations  19 
to  22.  In  these  two  cases,  respectively,  equations  19  and  29  for 
the  power  flow  yield 


W 


<Jl  = 


2ira‘ 


*1 


(97) 


where  ^  is  a  constant,  and,  noting  that  0^.  =  1  /J^~i 

W 


Ci  = 


/N 


(98) 


where  N  is  the  number  of  propagating  modes. 


In  contrast,  for  the  mode-block  approach,  the  physical 
interpretation  of  the  vector  components  and  transfer-matrix 
elements  provides  information  about  the  range  of  possible 
transmission  coefficients  for  a  component.  By  definition,  the 
components  of  the  mode  distribution  vector  u  are  the  amount  of 
power  flowing  in  a  mode,  so  u^>0  for  i  =  1,  .  ..,  N.  Since  the 
elements  of  the  scattering  matrix,  s,  represent  the  fraction  of 
power  scattered  from  one  mode  into  another,  0<s^j<l  for  i,j  =  1, 


given  by 

N 

W'  =  l  u'i  (99) 

i  =  1 


can  be  written  in  terms  of  the  input  mode  vector 

N 

W'  =  l  Ti  Ui  (100) 

i  =  1 


where  the  column  sum 


N 


Ti  =  .  I 


J  -  1 


li 


(101) 
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thus,  Tm±n  W  <  W'  <  W,  where  W  is  the  input  power  and  Tm^n 

and  ^raax  are  the  smallest  and  largest  column  sums  of  s. 
Consequently,  0  <  <1,  i  =  1,  ...»  N.  It  follows  from  the 

general  definition  of  a  reduce-dimension  transfer  matrix, 
equation  4,  that  the  column  sums 

M 

S-;  =  l  Hti  ,  j  =  1,  .  .  .  ,  M  (102) 

J  i  =  1  J 

can  be  rewritten  as 

N  M  +  N 

I  ,  M  T ik>  S)C1  n-j  -  l  7j  Tj j+  .  (103) 

k,  1  =  1  i  =  l  Jl  =  l 

Thus, 

^min  -  sj  -  ^max  (l°4) 

where  the  fact  that  the  column  sums  of  the  relevant  mode-block 

transformations  T  and  T+  are  all  equal  to  one  has  been  used. 

Since  0  <  H-^j ,  i,  j  =  1,  .  .  .  ,  M,  it  follows  from  equation  104  that 
0  s  j ,  Sj  <1,  which  is  consistent  with  the  interpretation  of  the 
elements  of  the  reduced-dimension  mode-block  transfer  matrix  as 
being  a  best  approximation  to  the  fraction  of  input  power  in  a 
mode-block  scattered  into  another  mode-block  on  output. 


Equation  100  can  be  rewritten  as 

W'  =  T( u)W  (105) 
where  TO u)  is  the  power  transmission  coefficient  of  the  component 
for  an  input  mode  distribution  u.  Thus, 


T{  u) 


N 

,  l  T i  Ui 

i  =  l 


(106) 
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Similarly  an  estimate  of  this  power  transmission  coefficient  is 

M 

7(c)  =  — -  (107) 

M 


where  by  equation  2  c  =  Tu.  Clearly,  Tm^n  <T  ( u)  ,  T  (c)  ^Tmax  and 
smin  *T(c)  <Smax,  where  Smin  and  Smax  are,  respectively,  the 
smallest  and  largest  column  sums  of  H. 

That  the  range  of  predicted  transmission  coefficients  can  be 
seen  immediately  from  the  column  sums  Sm^n  and  Smax  of  the 
reduced-dimension  transfer  matrix  in  mode-block  form  is  useful, 
but  this  does  not  provide  a  reliable  estimation  of  the  range  of 
observable  transmission  coefficients  determined  by  Tm±n  and  Tmax 
unless  the  partition  of  the  mode  distribution  has  been  made  into 
small  enough  blocks  so  that  it  is  unlikely  that  any  mode 
distribution  can  develop  with  power  predominantly  in  one  block, 
assuring  that  the  average  implied  in  equation  106  is  carried  out 
over  at  least  a  block.  This,  of  course,  makes  the  measurement  of 
a  mode  transfer  matrix  difficult  at  best. 

If  the  column  sums  Tm±n  =  7max,  then  also  the  column  sum  Sj  = 
S,  j  =  1,  . . . ,  M,  and  W'  =  S  W  and  such  an  element  would  be  an 

ideal  attenuator.  If  two  ideal  attenuators  are  concatenated,  the 
resulting  composite  system  will  also  be  an  ideal  attenuator  with 
S12  =  s1  S2.  It  would  appear  that  ideal  attenuators  can  be  fully 
characterized  by  a  scalar  transmission  coefficient  and  that 
systems  made  up  of  such  ideal  components  would  not  exhibit  the 
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concatenation  effect;  however,  this  would  be  true  if  only  power 
flow  were  of  interest.  If  the  signal  response  is  considered  as  in 
equation  48,  it  is  evident  that 

rmin  Smin  W  <  j  <  rmax  Smax  W  (108) 

where  rmax  and  rmin  are,  respectively,  the  largest  and  smallest 
photodetector-response-vector  components.  Therefore,  unless  the 
photodetector  is  characterized  by  a  uniform  vector  with  rm^n 
=  rmax,  the  observed  photocurrent  will  show  the  concatenation 
effect  even  though  it  may  be  true  that  Smin  =  Smax.  In  general, 
only  systems  made  up  of  components  whose  characterizing  transfer 
matrices  commute  will  not  be  subject  to  the  concatenation  effect. 

The  benefit  of  the  physical  interpretation  of  components 
represented  by  the  mode-block  formalism  can,  in  principle,  be 
extended  to  other  representations  through  a  transformation  as  in 
equations  87,  93,  and  94,  the  column  sums  then  being  evaluated  in 
the  resulting  mode-block  representation. 

EQUILIBRIUM  MODE  DISTRIBUTION 

The  eigenvectors,  vn,  and  eigenvalues,  \n,  of  scattering 
matrices  play  an  important  role  in  their  physical  interpretation. 
The  eigenvectors  of  a  scattering  matrix  of  a  component  represent 
the  mode-power  distributions  that  when  launched  into  a  component 
appear  unchanged  at  the  output  except  for  an  attenuation  factor  \ , 
the  associated  eigenvalue.  Consequently,  if  steady  state  or 
equilibrium  mode  distributions  exist  in  sequences  of  repeated 
components,  they  are  given  by  the  eigenvector  of  the  scattering 
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matrix  associated  with  the  largest  eigenvalue,  that  is,  greatest 
transmission  provided  it  represents  a  physically  realizable 
distribution. 

The  eigenvectors  and  eigenvalues  of  a  reduced-dimension 
transfer  matrix  should  have  a  similar  interpretation  to  those  of 
the  full  scattering  matrix.  However,  it  is  not  clear  that  enough 
information  is  retained  in  the  reduced  dimension  to  make  the 
predicted  equilibrium  distributions  useful.  That  is,  after  a 
number  of  passes,  say  k,  through  a  component  does  the  occurrence 
of  the  product  T+T  k-1  times,  see  equation  16,  diminish  the 
suitability  of  the  reduced-dimension  transfer  matrix  faster 
than  the  output  approaches  the  equilibrium  distribution?  The 
difficulty  here  is  analogous  to  those  that  arise  in  considerations 
of  the  concept  of  expressing  a  transfer  matrix  in  decibels,  which 
requires  that  the  matrix  be  written  as 

00 

H  =  exp(a)A  =  £  - -  An  (109) 

n  =  1  n! 

where  10  log  A  is  the  matrix  representation  of  H  in  decibels 
(references  12,  25,  and  29).  The  question  being,  does  the  error 
in  An  due  to  uncertainty  in  A  increase  more  rapidly  with  n  than 
the  series  converges?  Caution  should  be  exercised  when  using 
eigenvectors  of  reduced-dimension  transfer  matrices  to  predict 
equilibrium  mode  distributions. 


35 


COMPUTER  PROGRAMS 


Three  computer  programs  using  the  integral  form  of  the 
mode-power/near- field  relationship  for  converting 

near-field-pattern  data  into  reduced-dimension  transfer  matrices, 
detector  response  vectors,  or  source  vectors  have  been  written  in 
FORTRAN  and  are  available  on  IBM  compatible  MS-DOS  format  diskette 
from  the  author.  The  first  of  these  programs  uses  the  mode-block 
reduction.  The  other  two  reductions  are  based  on  orthogonal 
polynomial  expansions  of  the  mode-power  distribution.  One  of 
these  carries  out  the  expansion  to  a  large  number  of  terms,  say 
10,  and  then  truncates  the  expansion  at  the  number  of  terms  equal 
to  the  dimension  of  the  reduced-dimension  transfer  matrix,  while 
the  other  carries  out  the  expansion  only  to  the  desired  dimension 
of  the  reduced  matrix  from  the  start.  These  may  differ  somewhat 
in  the  final  result  primarily  because  the  weight  function  used  in 
the  least-square  fit  is  not  necessarily  the  same  as  that 
characterizing  the  orthogonal  polynomials,  but  also  because  the 
mode-power  distribution  appears  in  an  integral  in  the 
least-square-fit  error  function. 

In  all  three  programs,  alpha  index  of  refraction  profiles  and 
beta  mode  densities  were  assumed  (see  equations  32  and  33) ,  and 
the  choice  between  Chebyshev  and  Jacobi  orthogonal  polynomials 
with  weight  functions  l//x ( 1-x)  and  x,  respectively,  over  the 
domain  (0,1)  is  offered.  The  mode-block  program  uses  a 
mode-fraction  partition  so,  except  for  computation  error,  the 
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transfer  matrices  generated  by  it  should  be  independent  of  the 
choice  of  beta  (beta  must  not  be  zero)  as  well  as  the  selected  set 
of  polynomials.  The  other  two  programs  should  produce  results 
that  depend  on  these  choices.  In  all  three  cases,  the  value  used 
for  alpha  should  correspond  to  the  actual  index  of  refraction 
profile  of  the  fiber.  For  simplicity  the  programs  assume  that 
alpha  is  the  same  for  both  the  input  and  output  fibers.  However, 
the  core  diameter  as  well  as  the  normalized  frequency  may  be 
different  for  the  input  and  output  fibers.  Also,  source  and 
detector  vectors  are  referred  to  power-per-unit  area.  All  three 
programs  use  singular-value  decomposition  to  obtain  the 
pseudoinverse  of  the  measurement  matrix.  Consequently,  as  many 
measurement  sets  as  are  available  can  be  incorporated  into  a 
least-square  determination  of  a  transfer  matrix,  and  the  condition 
number  of  the  measurement  matrix  is  reported  to  permit  selection 
of  the  largest  acceptable  dimension  for  a  transfer  matrix. 

Also  available  on  diskette  from  the  author  are  two  programs  to 
determine  a  suitable  value  for  alpha  from  either  near-field 
pattern  data,  presuming  a  uniform  mode  distribution,  or  index  of 
refraction  profile  data.  In  the  first  case,  the  data  are  fitted 
by  a  function  of  the  form 

(1  -  p/a)a 

c1  -  +  c0  (110) 

71  -  K  (p/a)  2 

where  non-zero  K  accounts  for  leaky  modes  (reference  1),  and  in 
the  second  case,  by  a  function  of  the  form 

c2  [1  -  (p/a)a]  +  c1(p/a)a  +c0  .  (Ill) 
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The  coefficients  c 


0,  clf  and  c2  as  well  as  a  and  K  are  determined 
from  data  given  as  a  function  of  p/a. 

BANDWIDTH  CHARACTERIZATION 


The  transfer-matrix  formalism  presented  in  this  report  can  be 
extended  to  include  bandwidth  as  well  as  loss  characterization. 
Even  though  bandwidth  chracterization  has  not  yet  been  fully 
formulated  and  is  outside  the  scope  of  this  report,  a  brief 
description  of  a  possible  approach  for  doing  so  would  seem  in 
order  (reference  11) . 


The  basic  assumption  is  that  equation  1  can  be  extended  to 
include  a  time  dependence  in  the  mode-power  distribution  simply  by 
changing  the  form  to 


u.(t) 


00 

Sij  (t-T ) Uj ( r ) dr 

-00 


(112) 


where  now  0  <s^j(t)<l  is  to  be  interpreted  as  an  impulse  response 
for  scattering  of  power  from  the  jth  input  mode  into  the  ith 
output  mode. 


For  a  concatenated  system  consisting  of  two  elements 
characterized  by  s^1)  and  s(2) 

00 


00 

«:(t)-E  , 

1  r,j  =  1 


*00 


—00  —00 


Sir)(t-T2)  srj)(T2  ~  Ti)Uj (r1)dT1dr2  .(113) 
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By  taking  Fourier  transforms  of  both  sides  equation  113,  the  sums 
of  convolution  integrals  in  the  time  domain  may  be  replaced  by 
matrix  products  in  the  frequency  domain;  thus, 

u' (w)  =  s(2) (u)  u(o)  (114) 

where  the  Fourier  transform  relationship  is  indicated  by  replacing 
the  time  argument,  t,  by  the  frequency  argument,  o.  The  time 
dependent  quantities  u(t)  and  s(t)  are  real;  however,  in  general, 
the  frequency  dependent  Fourier  transforms  u(w)  and  s(w)  are 
complex  except  for  u  =  0,  where  they  are  real  and  correspond  to 
the  case  of  simple  loss  characterization  treated  in  the  bulk  of 
this  report. 


Just  as  in  the  case  of  loss  characterization,  a  frequency 
independent  transformation  T  can  be  introduced  to  define  a 
reduced-dimension  scattering  matrix  M(u>)  ,  where 

M(u)  =  T  s(w)  T+  .  (115) 

Also,  in  order  to  provide  a  measurable  approximation  H(w)  to  M(u) , 
a  frequency  dependent  continuous  mode  power  distribution  px(x,-j), 
which  can  be  related  to  the  Fourier  transform  of  the  time 


where  the  complex  coefficients  c^(w)  become  the  reduced  dimension 
representation  of  u(u>). 
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The  determination  of  the  complex  vector  c(w)  can  be  carried 
out  by  least-square  fitting  as  before  from  equation  51,  where  the 
matrix  A  is  real  and  identical  to  the  loss  characterization  case 
and  the  complex  vector  b(o)  must  be  found  from  measurements  of  the 
Fourier  transform  of  I(s,t);  thus 


fc>i(w)  = 


V‘ 


2na‘ 


s  I  ( s  ,0)  ) 


0^(x)dx  ds 


(118) 


f  (s) 


We  see  that  the  near-field  pattern  measurements  now  have  the 
additional  complication  of  requiring  phase  as  well  as  amplitude  to 
be  determined  at  each  point  in  the  pattern  over  the  expected 
modulation  bandwidth  of  the  system. 


Once  the  vector ization  process,  that  is,  the  determination  of 
b(w)  from  I(s,w),  is  complete,  then  input  and  output  data  matrices 
(c^  (w)  |  .  .  .  |  Cp(w)  )  and  (cj_ '  (w  )  |  .  .  .  |  cp  '  (w )  )  ,  respectively,  can  be 
formed  and  H(u)  determined  according  to  equation  63.  This 
determination  is  complicated  by  the  complex  vectors,  c(u), 
involved,  but  the  essential  question  is:  out  of  the  p  measured 
input  vectors  c^(o),  are  there  at  least  N,  where  N<P  is  the 
dimension  of  H,  independent  complex  vectors?  The  answer  is  yes  if 
an  ideal  chopped  modulation  is  assumed,  that  is,  it  is  sufficient 
that  the  input  near-field  patterns  can  be  written  in  the  separated 
form  I  (s,w)  =  I(s)F(<j),  where  F(w)  is  a  complex  function  of  unit 
amplitude.  The  M  vectors  b(w)  are  simply  real  vectors 
multiplied  by  a  complex  scalar,  and  the  resultant  c^(u>)  will  be 
independent  to  the  same  degree  as  could  be  achieved  in  simple  loss 
characterization. 
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Here  only  component  transfer  matrices  and  source  vectors  have 
been  mentioned,  but  the  detector  response  vector  r(w)  could  be 
similarly  determined  and  would,  in  general,  be  complex  as  would 
the  resulting  signal  current  j(w). 

Once  all  of  the  components  of  a  system  have  been 
characterized,  the  total  system  loss  and  bandwidth,  based  on 
consideration  of  |j(u)|2,  can  be  predicted  from  the  product  of 
component  matrices  and  the  detector-response  and  source  vectors. 

It  should  be  noted  that  even  matrices  that  have  only  real, 
frequency  independent  elements,  as  would  be  expected  of  small 
sized  passive  components,  must  be  included  in  the  product  as  they 
would  influence  not  only  the  loss  but  the  bandwidth  determined  in 
this  approach  as  well. 

If  a  mode-block  basis  is  used,  an  upper  limit  on  frequency 
response  of  a  system  can  be  expressed  in  terms  of  the  column  sums 
of  the  absolute  values  of  the  constituent  component  matrices. 
Consider,  for  example,  a  system  characterized  by 

j(w)  =  r  (w)  •  H(2)  (w)h(1)  (w)c(w)  .  (119) 

From  the  submultiplicative  property  of  the  maximum  column  sum 
matrix  norm  (reference  A-2,  section  5.6),  it  follows  that 

lj(«)l  *  |  |r(w)  |  |  x|  |H(2)  (W)  |  |  x|  IhU)  (w)  |  |  1\  |c(u)  |  |  1  (120) 

where,  for  an  M  by  M  matrix  A 

n 

I  I A |  | 1  =  max  l  |Aii |  .  (121) 

l<j<M  i  =  1  J 

Evidently  the  contribution  of  each  component  to  this  upper  limit 
can  be  obtained  independently;  however,  because  this  upper  limit 
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can  be  well  beyond  actual  performance,  caution  must  be  exercised 
in  using  such  limits. 


SUMMARY 

We  have  shown  that  a  reduced-dimension  transfer  matrix  can  be 
defined  in  terms  of  a  rectangular  matrix  transformation  from  the 
full-dimension  mode  space  to  a  reduced-dimension  space  and  its 
psuedoinverse .  Furthermore,  we  have  shown  that  the  reduced-dimension 
space  need  not  be  the  conventional  mode-block  partition.  Although 
the  reduced-dimension  transfer  matrices  are  in  some  sense  the  best 
possible  approximation  to  the  full  transfer  matrix,  as  a 
consequence  of  this  definition,  the  product  of  reduced  dimension 
transfer  matrices  inevitably  deviates  from  the  reduced  dimension 
representation  of  the  product  of  full  dimension  matrices. 

We  have  presented  a  least-square-fit  method  that,  to  a  large 
extent,  circumvents  errors  resulting  from  having  to  compute  the 
derivative  of  noisy  data,  as  is  required  in  the  conventional 
method,  and  permits  the  natural  inclusion  of  additional  data  sets 
beyond  the  minimum  number  of  the  dimension  of  the  reduced 
representation.  As  a  consequence  of  using  singular-value 
decomposition  in  this  method,  the  condition  number  of  the  data 
matrix  is  obtained  and  provides  a  measure  of  the  independence  of 
the  launch  conditions,  which  is  essential  to  the  success  of 
measurements  of  transfer  matrices. 
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We  have  shown  that  all  reduced-dimension  representations  are 
equivalent  in  the  sense  that  one  can  be  obtained  from  another 
without  knowledge  of  the  full  representation,  although  to  do  so 
results  in  some  unavoidable  loss  in  information.  Consequently, 
valid  reduced-dimension  representations  should  be  attainable 
directly  from  vectorization  of  the  near-field  patterns  without 
first  developing  the  mode  distribution. 

In  addition  to  the  consideration  of  loss  characterization  by 
transfer  matrix,  we  have  briefly  touched  upon  the  extension  of  the 
transfer  matrix  approach  to  bandwidth  characterization.  This 
extension  needs  both  experimental  and  theoretical  development  to  be 
useful  for  applications. 

Computer  programs  written  in  FORTRAN-77  for  the  computation  of 
various  representations  of  reduced  dimension  transfer  matrices 
from  near-field  patterns  are  available  from  the  author. 
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APPENDIX  A.  MATHEMATICAL  BACKGROUND 

GENERALIZED  INVERSE  (references  A1  and  A2) 

The  least-square  problem  to  find  an  x  to  minimize  | |Ax-b| I2 / 
where  A  is  an  m  by  n  matrix,  x  and  b  are  n  and  m  component 
vectors,  respectively,  and  ||«||2  is  the  Euclidean  vector  norm, 
has  the  unique  solution 

xLS  =  ( AhA) _1AHb  (Al) 

where  A^  is  the  Hermitian  transpose  of  A,  if  the  rank  of  A  is  n, 
that  is,  if  the  columns  of  A  are  linearly  independent.  Otherwise, 
the  vector  xLS  such  that  | |x| |2  is  minimized  among  all  those  x 
that  produce  the  minimal  | | Ax— b { | 2  may  be  obtained  from 


XLS  =  A+b  (A2) 

where  A+,  an  n  by  m  matrix,  is  typically  defined  to  be  the  unique 

matrix  that  satisfied  the  Moore-Penrose  conditions 

A  A+  A  =  A  (A3) 

A+  A  A+  =  A+  ( A4 ) 

(A  A+)H  =  A  A+  (A5) 

(A+  A)H  =  A+  A  ( A6 ) 


and  is  known  as  the  Moore-Penrose  generalized  or  pseudoinverse  of 
A.  A  property  of  particular  importance  here  is  that  X  =  A+  is  the 
unique  matrix  that  minimizes  the  Frobinius  matrix  norm  |  I  AX— Im |  |p, 
where  Im  is  the  m  by  m  unit  matrix. 
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SINGULAR-VALUE  DECOMPOSITION  (references  A3  and  A4) 


Any  m  by  n  matrix  A  of  rank  r  can  be  decomposed  into  the 
factored  form 

A  =  U  X  VH  (A7  ) 

where  U  and  V  are  m  by  m  and  n  by  n  unitary  matrices  respectively 
and  ][  is  an  m  by  n  diagonal  matrix  with  elements  or  such  that 

a1>a2a>. . ->ar+1  =  ...  =  aq  =  0  (A8) 

in  which  q  is  the  lesser  of  m  and  n.  The  ar  are  uniquely 
determined  by  A,  being  the  square  root  of  the  eigenvalues  of  AAH, 
and  the  non-zero  ar  are  known  as  the  singular  values  of  A.  In 
terms  of  the  singular-value  decomposition,  the  generalized  inverse 
of  A  can  be  written 

A+  =  V  l+  UH  ( A9 ) 

where  £+  is  a  n  by  m  diagonal  matrix  with  elements  d^  such  that 

1 

<*i  =  ~  ,  i  =  l,  ...,r  .  (A10) 

ai 

Thus,  the  singular-value  decomposition  provides  a  means  to 
determine  A+.  Stable  algorithms  for  carrying  out  least-square 
solutions  by  singular-value  decomposition  without  explicit 
determination  of  the  eigenvalues  of  AAH  are  available  and  only 
require  about  three  times  the  computational  effort  as  other 
least-square  approaches  (reference  A5) . 

LEAST-SQUARE  PERTURBATION  (references  Al  and  A5) 

Suppose  we  have  solutions  Xj-^  and  yLS  to  the  least-square 
problems  |  |  Ax-b  |  \2  minimized  and  |  |  ( A+<5A)  y-  (b+<5b)  |  \2  minimized, 
where  both  A  and  its  perturbation  A+<5A  are  m  by  n  matrices  with 


A-2 


m>n  and  of  full  rank  n  and  5  and  its  perturbation  <5b  are  m 
dimensional  vectors  with  b  0.  The  difference  in  these  two 
solutions  satisfies  to  first  order  in  £  the  relationship 

I  I  yLS~xLS  I  I  2  .  ..  .  -  . 

-  <  £  K  (2  +  K  sm  e)/cos  ©  (All) 

I  I  XLS  I  I  2 


where  we  assume  that 


£ 


max 


~IUa||2 

1  I  All  2 


and  that 

I  I AxLS"bl  I  2 

sin  ©  =  -  *  1 

Mb||2 


(A12) 


(A13) 


with  K,  the  condition  number  of  A,  equal  to  |  | A |  |2  I  I A+ 1  | 2  =  °i/°n 
the  ratio  of  the  maximum  to  minimum  singular  values  of  A. 

Equation  All  gives  an  upper  bound  on  the  difference  between  the 
solutions  of  a  least-square  problem  and  its  perturbation  and  as 
such  may  be  much  larger  than  the  actual  difference.  However,  it 
would  be  wise  to  be  cautious  when  large  values  of  condition  number 
occur  because,  in  contrast  to  the  linear  equation  case,  where  m=n, 
differences  can  be  magnified  by  the  square  of  K  for  poorly  fitting 
solutions.  Furthermore,  the  magnification  could  be  even  greater 
than  indicated  by  equation  All  for  situations  where  A  is  not  of 
full  rank. 
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